function pmpv = pmp_Sample(hfile, coords, mnum)
% PMP::Sample  - sample PMP at given coordinates
%
% FORMAT:       pmpv = pmp.Sample(coords [, mnum]);
%
% Input fields:
%
%       coords      Cx3 cartesian coordinates to sample PMP at (centered!)
%       mnum        optional map number (default: 1)
%
% Output fields:
%
%       pmpv        Cx1 PMP values at coords

% Version:  v0.7f
% Build:    9030320
% Date:     Mar-03 2009, 8:14 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'pmp') || ...
   ~isa(coords, 'double') || ...
    ndims(coords) > 2 || ...
    size(coords, 2) ~= 3
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get file contents
bc = bvqxfile_getcont(hfile.L);
r1 = bc.ThetaResolution;
r2 = bc.PhiResolution;
if nargin < 3 || ...
   ~isa(mnum, 'double') || ...
    numel(mnum) ~= 1 || ...
    isinf(mnum) || ...
    isnan(mnum) || ...
    mnum < 1 || ...
    mnum > numel(bc.Map)
    mnum = 1;
else
    mnum = round(mnum);
end

% get spherical coordinates
scoords = spherecoords(coords);

% compute indices
crd1 = 1 + mod(r1 + (r1 / (2 * pi)) * scoords(:, 3), r1);
crd2 = 1 + mod(r2 + (r2 / pi) * scoords(:, 2), r2);

% build map for interpolation
imap = bc.Map(mnum).PMPData([1:r1, 1], [1:r2, r2]);

% interpolate
pmpv = flexinterpn_method(imap, [crd1, crd2], 0, 'linear');
